#import rasterio
import numpy as np
import scipy.io as sio
import h5py
import gzip
import glob
import os
import time

#Begin timer.
start = time.time()
print 'START'

#PATHS
#dirbl: Directory where BIL format files downloaded from PRISM server are stored.
dirbl='/Volumes/TOSHIBA/climgroup_research_macbookair/HWA/manuscripts/C__HWAmicroclimate/companion_dataset/05_prism/bil'
#dirble: Directory where PRISM normals files are stored.  PRISM elevation variable is stored in PRISM normals file "PRISM_us_dem_4km_bil.bil". 
dirble='/Volumes/TOSHIBA/climgroup_research_macbookair/HWA/manuscripts/C__HWAmicroclimate/companion_dataset/05_prism/1991-2020_normals/bil' 
#dirh5: Directory where HDF-5 format files are to be stored.
dirh5='/Volumes/TOSHIBA/climgroup_research_macbookair/HWA/manuscripts/C__HWAmicroclimate/companion_dataset/05_prism/h5/'

#SUBSET FULL PRISM CONUS DOMAIN:
ss_w=-98.6;
ss_e=-71.3;
ss_s=36.9;
ss_n=49.6;

#####STEP 1: STATIC VARIABLES######
print 'STEP 1: STATIC VARIABLES'
#CREATE STATIC ARRAYS AND WRITE TO HDF-5 FILE
#(1): LAT,LON
#Construct latitude and longitude arrays
#Contents of *bil.hdr file.
BYTEORDER='I';
LAYOUT='BIL';
NROWS=621;
NCOLS=1405;
NBANDS=1;
NBITS=32;
BANDROWBYTES=5620;
TOTALROWBYTES=5620;
PIXELTYPE='FLOAT';
ULXMAP=-125;
ULYMAP=49.9166666666687;
XDIM=0.0416666666667;
YDIM=0.0416666666667;
NODATA=-9999;
#From : http://desktop.arcgis.com/en/arcmap/10.3/manage-data/raster-and-images/bil-bip-and-bsq-raster-files.htm
#ulxmap:The x-axis map coordinate of the center of the upper left pixel. If you specify this parameter, set ulymap, too; otherwise, a default value is used.
#ulymap:The y-axis map coordinate of the center of the upper left pixel. If this parameter is specified, ulxmap must also be set; otherwise, a default value is used.
temlats= np.empty((NROWS,NCOLS), dtype=np.float, order='F');temlats[:]=np.nan;
temlons= np.empty((NROWS,NCOLS), dtype=np.float, order='F');temlons[:]=np.nan;
for j in range(NROWS-1,-1,-1):
    temlats[j,:]=ULYMAP+((j-NROWS+1)*YDIM);#7/31/19: Replaced NROWS with NROWS+1
for i in range(0,NCOLS,1):
    temlons[:,i]=ULXMAP+((i)*XDIM);#7/31/19: Replaced i-1 with i
I=(temlats>=ss_s)&(temlats<=ss_n)&(temlons>=ss_w)&(temlons<=ss_e);
jll=np.where(I>0)[0][0]
ill=np.where(I>0)[-1][0]
jur=np.where(I>0)[0][-1]
iur=np.where(I>0)[-1][-1]
LATS=temlats[jll:jur+1,ill:iur+1];LONS=temlons[jll:jur+1,ill:iur+1]
###(2): MASK
###Apply HWA state mask (only show values for MN,WI,MI,IL,IN,OH,PA,NY)
###See "make_mask.m" for details of mask construction.
###MASK = sio.loadmat(dirmt+'/hwastatemask.mat')['mask'];
#(3): ELEV
filename=dirble+'/elev/PRISM_us_dem_4km_bil.bil'
prism_path = filename
tem_array = np.fromfile(prism_path, dtype=np.intc)
tem_array = tem_array.reshape(NROWS, NCOLS)
tem_array[tem_array == NODATA] = 99999
tem_array=np.flipud(tem_array)
ELEV=tem_array[jll:jur+1,ill:iur+1]
#WRITE ALL STATIC VARIABLES TO HDF-5 FILE
print 'Writing ELEV,LATS,LONS to static HDF-5 file'
with h5py.File(dirh5+'/PRISM_D2_static.h5', 'w') as hf:
	hf.create_dataset('ELEV',data=ELEV,compression="gzip", compression_opts=9)	
	hf.create_dataset('LATS',data=LATS,compression="gzip", compression_opts=9)		
	hf.create_dataset('LONS',data=LONS,compression="gzip", compression_opts=9)	
	
#####STEP 2: TIME-DEPENDENT VARIABLES######
print 'STEP 2: TIME-DEPENDENT VARIABLES'
for yr in range(1981,2021,1):
	print 'Processing '+str(yr)	
	#TMIN
	numfiles=np.size(sorted(glob.iglob(dirbl+'/tmin/PRISM_tmin_stable_4kmD2_'+str(yr)+'*')))
	TMIN= np.empty((np.size(LATS,0),np.size(LATS,1),numfiles), dtype=np.float, order='F');TMIN[:]=np.nan;	
	d=0
	for filename in sorted(glob.iglob(dirbl+'/tmin/PRISM_tmin_stable_4kmD2_'+str(yr)+'*')):
		prism_path = filename
		tem_array = np.fromfile(prism_path, dtype=np.float32)
		tem_array = tem_array.reshape(NROWS, NCOLS)
		tem_array[tem_array == NODATA] = np.nan
		tem_array=np.flipud(tem_array)
		TMIN[:,:,d]=tem_array[jll:jur+1,ill:iur+1]
		d=d+1
	print 'Writing TMIN to single-year HDF-5 file'		
	with h5py.File(dirh5+'/PRISM_D2_'+str(yr)+'.h5', 'a') as hf:
		hf.create_dataset('TMIN',data=TMIN,compression="gzip", compression_opts=9)	
	del TMIN							
print 'ALL DONE!'
#Stop timer.
end = time.time()
print(end - start)
